package scu.maqiang.cte;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class BrickResultT4_SC2 {

	public static void main(String[] args) {
	    double kxx = 1.0, density = 1.0, specficHeat = 1.0;
	    double E = 2.0, Nu = 0.3, alpha = 0.05;
	    double beta = E / (1 - 2 * Nu)  * alpha;

	    double theta = 0.5, endTime = 10, dt = endTime / 1000;
	    double omega = 0.25;

	    Mesh3T4 mesh3D = new Mesh3T4().cube(10, 10, 20);
	    mesh3D.scale(0.5, 0.5, 1.0);

	    //�����¶ȳ���Ԫ�ռ�, ��װ��Ӧ�������
	    FES3T41 fs = new FES3T41(mesh3D);
	    SRMatrix LMatrix = new SRMatrix(fs.GetNdof());
	    SRMatrix RMatrix = new SRMatrix(fs.GetNdof());
	    ScalarFunc kxxFunc = (xy, llabel, param) -> kxx;
	    ScalarFunc rhoFunc = (xy, llabel, param) -> density * specficHeat;
	    fs.assembleHeatThetaScheme(new ScalarFunc[]{kxxFunc, rhoFunc}, 
	    		                   new double[] {dt, theta}, 
	    		                   new BVPType[] {BVPType.COMMON, BVPType.CONSIST_MASS},
	    		                   LMatrix, RMatrix);
	    double[] x = new double[fs.GetNdof()];
	    double[] RHS = new double[fs.GetNdof()];

	    //���õ��Գ�����Ԫ�ռ�, ��װ��Ӧ�������
	    FES3T42 fs2 = new FES3T42(mesh3D);
	    SRMatrix K = new SRMatrix(fs2.GetNdof());
	    SRMatrix M = new SRMatrix(fs2.GetNdof());
	    double[] d = new double[fs2.GetNdof()];
	    double[] d0 = new double[fs2.GetNdof()];
	    double[] vel = new double[fs2.GetNdof()];
	    double[] vel0 = new double[fs2.GetNdof()];
	    double[] acc = new double[fs2.GetNdof()];
	    double[] acc0 = new double[fs2.GetNdof()];
	    double[] heatLoad = new double[fs2.GetNdof()];
	    double[] totalLoad = new double[fs2.GetNdof()];
	    
	    fs2.assembleStiff(new double[] {E, Nu}, BVPType.COMMON, K);
	    fs2.assembleMass(new double[] {density}, BVPType.CONSIST_MASS, M);
	    ScalarFunc myFunc = (coord, llabel, param) -> coord[2] > 0.99999? 1.0 : 0.0;
	    double[] x0 = fs.valueFromFunc(myFunc, null);
	    ScalarFunc betaFunc = (xy, llabel, param) -> beta;
	    
	    double[][] vecX = new double[1][];
	    vecX[0] = x0;
	    fs2.assembleHeatLoad(new ScalarFunc[] {betaFunc}, null, vecX, BVPType.COMMON, heatLoad);
	    
	    //求解初始加速度
	    NewIterSSolver solverM = new NewIterSSolver(M);
	    solverM.PCGSSOR(heatLoad, acc, 1.5, 1);
	   
	    for(int i = 0, n = fs2.GetNdof(); i < n; i++) {
	    	acc0[i] = acc[i];
	    }

	    double[][] uvw = new double[3][fs.GetNdof()];
	    fs2.extractComponent(acc0, uvw);
	    
	    String dir = "BrickResultT4\\";

	    mesh3D.toTecplot(dir + "acc.dat", uvw);
	    MVO.fill(uvw, 0.0);
	    
	    //K.sort();
	    //System.out.println(K);
	    //System.exit(0);
	    
	    SRMatrix ELMatrix = K.add((1.0 / (0.25 * dt * dt)), M);
	    //ELMatrix.sort();
	    //ELMatrix.toString();
	    //System.out.println(ELMatrix);
	    //System.exit(0);
	    
	    double time = 0.0;
	    
	    String heatoutputFile = dir + "Heat100_2.dat";
	    mesh3D.toTecplot(heatoutputFile, time, x);

	    String elasticOutputFile =  dir + "Deform100_2.dat";
	    mesh3D.toTecplotDeform(elasticOutputFile, 0.0, uvw, 0.0);
	    //call mesh3D.exportMeshDeformTime(elasticOutputFile, time, x, x, 0.0D0)
	    
	    int lineNum = 100;
	    double dis = 1.0 / lineNum;
	    double[] lineX = new double[lineNum + 1];
	    double[] lineY = new double[lineNum + 1];
	    Arrays.setAll(lineX, i -> i * dis);

	    String XYLineFile = dir + "Line100_2.dat";
	    Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);

	    //x = 1.0D0
	    String expansion = dir + "Expansion100_2.dat";
	    String temp00File = dir + "temp00_100_2.dat";
	    FileWriter fw = null;
	    FileWriter fw2 = null;
	    try {
			fw = new FileWriter(expansion);
			fw.write("0.0\t0.0\n");
			fw2 = new FileWriter(temp00File);
			fw2.write("0.0\t0.0\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	    //XYStressYFile = trim(dir)//'StressYLine100_2.dat'
	    //call exportXYTime(XYStressYFile, lineX, time, lineY)
	    
	    fs.applyBC_MBN(LMatrix, Direct.All, 2);
	    
	    fs2.applyBC_MBN(K, Direct.X, 3, 4, 5, 6);
	    fs2.applyBC_MBN(K, Direct.Y, 3, 4, 5, 6);
	    fs2.applyBC_MBN(K, Direct.Z, 1);
	    
	    fs2.applyBC_MBN(ELMatrix, Direct.X, 3, 4, 5, 6);
	    fs2.applyBC_MBN(ELMatrix, Direct.Y, 3, 4, 5, 6);
	    fs2.applyBC_MBN(ELMatrix, Direct.Z, 1);
	    

	    
	    NewIterSSolver SolverLM = new NewIterSSolver(LMatrix);
	    double[] Md = new double[fs2.GetNdof()];
	    double[] MMd = Md.clone();
	    
	    NewIterSSolver solverEL = new NewIterSSolver(ELMatrix);
	    time = 0.0;
	    double coef = (1 - Nu) / (alpha * (1 + Nu));
	    while(time < endTime) {
	        time += dt;
	        System.out.println("time: " + time);
	        
	        fs.applyBC_MBN(RHS, Direct.All, 1.0, 2);
	        SolverLM.PCGSSOR(RHS, x, 1.5, 1);
	        mesh3D.toTecplot(heatoutputFile, time, x);
	        //RHS = RMatrix * x
	        RMatrix.mul(x, RHS);
	        
	        for(int i = 0; i < lineNum + 1; i++) {
	        	lineY[i] = fs.getValue(x, new double[] {0.25, 0.25, lineX[i]});
	        }
	        Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);
	        

	        Arrays.fill(heatLoad, 0.0);
	        vecX[0] = x;
	        fs2.assembleHeatLoad(new ScalarFunc[] {betaFunc}, null, vecX, BVPType.COMMON, heatLoad);
	        //call fs2.assembleHeatLoad(beta, x, 0, HeatLoad)
	        //M.mul(1.0 / (0.25 * dt * dt), d, Md);
	        Arrays.setAll(Md, i -> 1.0 / (0.25 * dt * dt) * d[i] + (1.0 / (0.25 * dt)) * vel[i] + (0.5 / 0.25 - 1.0) * acc[i]);
	        M.mul(Md, MMd);
	        for(int ii = 0, nn = fs2.GetNdof(); ii < nn; ii++) {
	        	totalLoad[ii] = heatLoad[ii] + MMd[ii];
	        }
	        //System.out.println(MVO.toString(totalLoad));
	        //System.out.println(MVO.max(totalLoad) + "\t" + MVO.min(totalLoad));
	        
		    fs2.applyBCZero(totalLoad, Direct.X, 3, 4, 5, 6);
		    fs2.applyBCZero(totalLoad, Direct.Y, 3, 4, 5, 6);
		    fs2.applyBCZero(totalLoad, Direct.Z, 1);

		    solverEL.PCGSSOR(totalLoad, d, 1.5, 1);
		    fs2.extractComponent(d, uvw);
		    mesh3D.toTecplotDeform(elasticOutputFile, time, uvw, coef);
		    
	        //call fs2.computeStrain(u, v, ex, ey, exy)
	        //call fs2.computeStress(ex, ey, exy, u, E, Nu, tx, ty, txy, u, 1)

	        //do i = 1, lineNum + 1
	        //    call fs.valueXY(ty, 0.25D0, lineX(i), val)
	        //    lineY(i) = val 
	        //end do
	        //call exportXYTime(XYStressYFile, lineX, time, lineY)
		    
		    Arrays.setAll(acc, i -> 1/(omega*dt*dt) * (d[i] - d0[i]) -1/(omega * dt) * vel0[i] - (0.5 /omega - 1) * acc0[i]);
		    Arrays.setAll(vel, i -> vel0[i] + ((1-0.5) * acc0[i] + 0.5 * acc[i]) * dt);
		    Arrays.setAll(acc0, i -> acc[i]);
		    Arrays.setAll(vel0, i -> vel[i]);
		    Arrays.setAll(d0, i -> d[i]);
	        
		    try {
				fw.write(time + "\t " + fs.getValue(uvw[2], new double[] {0.25, 0.25, 1.0}) * coef + "\n");
				fw2.write(time + "\t " + fs.getValue(x, new double[] {0.0, 0.0, 0.0}) + "\n");
			} catch (Exception e) {
				e.printStackTrace();
			}
	    }
	    try {
	    	fw.flush();
			fw.close();
			fw2.flush();
			fw2.close();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	}

}
